Microcompartment assembly around multicomponent fluid cargoes

This article describes dynamical simulations of the assembly of an icosahedral protein shell around a bicomponent fluid cargo. Our simulations are motivated by bacterial microcompartments, which are protein shells found in bacteria that assemble around a complex of enzymes and other components involved in certain metabolic processes. The simulations demonstrate that the relative interaction strengths among the different cargo species play a key role in determining the amount of each species that is encapsulated, their spatial organization, and the nature of the shell assembly pathways. However, the shell protein–shell protein and shell protein–cargo component interactions that help drive assembly and encapsulation also influence cargo composition within certain parameter regimes. These behaviors are governed by a combination of thermodynamic and kinetic effects. In addition to elucidating how natural microcompartments encapsulate multiple components involved within reaction cascades, these results have implications for efforts in synthetic biology to colocalize alternative sets of molecules within microcompartments to accelerate specific reactions. More broadly, the results suggest that coupling between self-assembly and multicomponent liquid–liquid phase separation may play a role in the organization of the cellular cytoplasm.


I. INTRODUCTION
Compartmentalization is essential for many biological functions, including metabolism, cellular signaling, and genetic storage. While membrane-enveloped organelles are a prominent mode of compartmentalization in eukaryotic cells, it has become apparent that liquid-liquid phase separation [1][2][3][4][5][6][7][8][9][10][11][12][13][14] and proteinaceous organelles [15][16][17] play an important role in organizing the cytoplasm for all kingdoms of life. For example, in bacteria, certain metabolic pathways are enabled by bacterial microcompartments, which are organelles consisting of a large protein shell assembled around a dense complex of enzymes and reactants. [18][19][20][21][22][23][24][25] Other protein-shelled compartments include encapsulins 26,27 and gas vesicles 26,28 in bacteria and archaea and vault particles in eukaryotes. 29 From a practical perspective, there is great interest in exploiting principles and materials from biology to achieve similar nanoscale compositional control for synthetic biology and drug delivery applications. In particular, researchers have demonstrated the ability to target new molecules to microcompartment interiors and to transfect the systems into non-native organisms, including bacterial and plant cells, suggesting a basis for designing microcompartments as customizable nanoreactors (e.g., .
Although previous experiments successfully targeted particular cargo species to microcompartment interiors, 22,[30][31][32][33][34][35][36][37][38][39][40][41][42][43] the factors that control cargo coalescence and encapsulation remain incompletely understood. In some microcompartments, shell-cargo attractions are mediated by "scaffold proteins" such as the CcmN protein in β-carboxysomes, 120 while in other systems core enzymes have a short "encapsulation peptide" sequence that binds to the shell inner surface. 26,31,[35][36][37]39,121 Similarly, cargo-cargo attractions may be mediated by scaffolds, such as CcmM in β-carboxysomes 120 or potentially through direct pair interactions between cargo molecules. 36,37 Previous experimental and modeling studies suggest that the packaged cargo of microcompartments undergoes phase separation, either prior to or during assembly of the outer protein shell. [60][61][62]65,120,[122][123][124][125][126] Furthermore, some microcompartments fulfill similar functions as liquid-liquid phase separated domains in cells-for example, the composition, structure, and function of the RuBisCO complex within carboxysomes have strong similarities to those of the pyrenoid, which is a liquid domain consisting of RuBisCO and other components that enables carbon fixation within plant cells. [127][128][129][130] It is, therefore, likely that some factors that control the composition of liquid-liquid phase separated domains in cells (e.g., Refs. 2, 6, 8, 10, 12-14, and 131-145) also affect the composition of microcompartment cargoes. However, in the case of microcompartment assembly, thermodynamic and kinetic factors resulting from coupling between liquid-liquid phase separation and assembly of the crystalline shell could have significant effects. Previous computational modeling of microcompartment assembly has focused on encapsulation of a single cargo species, either driven by direct cargo-cargo pair attractions 60,61,65 or mediated by scaffolds. 62 In this work, we build on these previous studies, by studying a minimal computational model for a microcompartment shell that assembles around a cargo containing two species. We investigate the amount and composition of the encapsulated cargo as a function of the affinities among these cargo species, and of the shell-cargo and shell-shell affinities. The simulation results suggest that the relative cargo-cargo affinities are the most important factor in determining the amount and composition of packaged cargo. However, within certain parameter ranges, varying the shell-cargo and shell-shell affinities provides an alternative means to sensitively control the packaged cargo. These effects are governed by a combination of thermodynamic and kinetic factors arising due to coupling between cargo coalescence and shell assembly.

Interaction parameters
Efficient assembly and cargo loading generically require three classes of interactions: interactions between shell subunits that drive shell assembly, interactions between shell subunits and one or more cargo species, and interactions between cargo particles that drive cargo coalescence. 60 While in some systems, one or more of these interactions can be mediated by auxiliary proteins (e.g., the CCmM and CCmN proteins in β-carboxysomes), for simplicity we model all interactions as direct pair interactions in this work.
We have explored assembly dynamics over a range of the simulation control parameters that most strongly affect the amount and composition of encapsulated cargo-these are the interaction potential well-depths that control the binding free energy (affinity) between each type of cargo pair, ε GG , ε RR , and ε RG , as well as between shell-cargo and shell-shell pairs, ε SC and ε SS (see Appendix B). To minimize the number of parameters, we assume that both cargo species have the same binding affinity to shell subunits; i.e., ε SC is equal for R and G particles. For convenience, throughout this article, we refer to the well-depths as "affinities." We approximately relate ε SC and ε SS to the actual affinities in Appendix A and Table I. Throughout this article, all energy values are given in units of the thermal energy, k B T, and all lengths are given in units of the cargo diameter, r * , which translates to r * ≈ 13 nm for carboxysomes (see Sec. IV B). We consider ranges of affinity parameters relevant to biological microcompartments (see Sec. IV B for details). Except where mentioned otherwise, we consider two shell-shell affinities (ε SS /k B T = 2, 3.5), which are, respectively, below and above the transition between one-step and two-step assembly pathways (discussed in Sec. II C). Both cases are below the threshold for assembly of empty shells (containing no cargo). Importantly, dimerization is very unfavorable, but the assembled shells are stable (see Table I). The unfavorable dimerization affinities reflect the essential role that cargo plays in assembly at these shell-shell affinities, the cooperative nature of assembly, and the fact that a significant nucleation barrier is essential to avoid kinetic traps. 60,80,89,157 We primarily focus on a representative shell-cargo affinity ε SC = 6k B T, which corresponds to a free energy per shell subunit of g SC ≈ −14k B T. Our simulations span the range of cargo-cargo affinities in the vicinity of vapor-liquid coexistence (ε vl ≈ 1.3k B T at our simulated cargo concentration), but they remain below the threshold for cargo crystallization at ε fs ≈ 3k B T. With ε vl as a reference point, we denote cargo-cargo affinities as "weak" for ε RR ≲ 1.1k B T, "moderate" for 1.2 ≲ ε RR ≲ 1.4, and "strong" for ε RR ≳ 1.5k B T (with analogous definitions for ε GG and ε RG ).
Our shell subunit and cargo concentrations map to ∼1 and 10 μM, respectively, which are reasonable based on quantities found within bacterial cells.

Simulated systems
There are multiple thermodynamic and kinetic factors that could influence the amount and composition of encapsulated cargo. Thermodynamic effects include the interplay between cargo-cargo and shell-cargo affinities, and the finite size of the encapsulated cargo globule. Kinetic effects arise from the fact that once a shell closes, even with relatively weak shell-shell affinities typical of productive assembly, exchange of shell subunits or cargo particles does not occur on experimentally relevant timescales. Thus, an assembling shell may trap its contents far out of equilibrium. To distinguish these effects, we begin by describing large simulations of cargo in the absence of shell subunits to approximate a bulk system (Sec. II B). We then describe the dynamical assembly of shells and cargo particles (Secs. II C and II D). We compare the results of these finite-time dynamical simulations against simulations with assembled but permeable shells, which allow the encapsulated cargo to equilibrate with the bulk.

B. Cargo phase behavior without shell subunits
We begin by briefly summarizing the bulk phase behavior of the cargo in the absence of shell subunits, as a reference point from which to understand how shell assembly can change cargo coalescence. We focus on parameter regimes relevant to our shell assembly simulations; more comprehensive descriptions of the phase behavior of a binary Lennard-Jones (LJ) system can be found in, e.g., Refs. 158-160. In particular, we consider equal stoichiometry between R and G molecules and cargo-cargo affinities (εij < 1.8 for ij = RR, GG, RG) that maintain the system in vapor and/or liquid phases but are below the crystallization threshold (the vapor-fluid and crystallization transitions for a single species occur at ε vl = 1.3 and ε fs ≈ 3). Within these limits, the phase behavior can be classified as follows: homogeneous demixed (no phase coexistence), phase separation of one cargo species (a dense phase rich in one cargo species in coexistence with a dilute phase containing an excess of the other), phase separation and mixing of both cargo species (coexistence between dilute and dense phases, with both species homogeneously mixed within both phases), and phase separation with demixing (coexistence between three phases-a dilute phase, a dense phase rich in R, and a dense phase rich in G). We consider weak, moderate, and strong R-G affinities: ε RG = 1.0, 1.3, 1.6. For each of these, Fig. 2 summarizes the bulk phase behaviors as a function of the R-R and G-G affinities ε RR and ε GG . Specifically, Figs. 2(a)-2(c) compare the driving force for phase separation for each of the two species by showing the fraction of R particles f R in the high-density phase (hereafter referred to as "globule"). Figures 2(d)-2(f) quantify the extent of nonrandom compositional mixing within the high-density phase, f UL . Here, we define the mean fraction of unlike neighbors in the first solvation shell around each particle asf UL ≡ n RG n RR +n GG n RG , with nij the number of cargo neighbor pairs of species i and j with a neighbor pair defined as two particles with separations of r < 1.2. For random mixing, we expect f rand = 2 f R (1 − f R ), so the plots indicate the difference between randomness and the results: For weak R-G affinities [ε RG = 1.0, Fig. 2(a)], the system exhibits no phase separation, coexistence between a phase rich in one species and a dilute phase, or coexistence among a dilute phase and separate phases for each species. In the latter case, the R-rich and G-rich phases are attracted to each other by the weak but nonzero R-G affinities, leading to formation of one globule with an interface separating the two phases. The interface between these two phases becomes more diffuse with decreasing ε RR and ε GG and/or increasing ε RG because the surface tension between the two domains decreases. For moderate R-G affinities, i.e., ε RG = 1.3, Fig. 2(b), we observe no phase separation, coexistence between a phase rich in one species and a dilute phase, or coexistence between a dilute phase and a dense phase containing both species.
For strong R-G affinities [ε RG = 1.6, Fig. 2(c)], we always observe mixing of both species within the dense phase, with the fraction of R or G particles depending on the relative magnitudes of ε RR and ε GG . Since the primary driving force for phase separation is ε RG in this case, fluctuations in f R and f UL are smaller than for the other systems.

C. Assembly pathways
We now consider how cargo and shell affinities affect the dynamics of assembly and cargo encapsulation. In a previous work, simulations with a single cargo species showed that assembly pathways are strongly affected by the net strength of cargo-cargo affinities. 60 In particular, at our simulated volume fraction of 0.003, weak cargo-cargo affinities (ε RR < 1.7) and/or the higher shell-shell affinities (ε SS > 2.5) lead to one-step pathways in which shell assembly is concomitant with cargo coalescence, while strong cargo-cargo and lower shell-shell affinities (ε RR > 1.7 and ε SS < 2.5) allow two-step pathways in which a cargo globule coalesces, followed by disordered adsorption of shell subunits onto the globule surface and then cooperative rearrangement into an assembled shell. Other interaction combinations, such as high cargo-cargo and shell-shell affinities, lead to pathways in between these two extremes.
The case of two cargo species exhibits qualitatively similar classes of assembly pathways. Figures 3(a) and 3(b) show snapshots from typical simulations at shell-shell affinity ε SS = 3.5, leading to one-step assembly pathways. In Fig. 3(a), the cargo-cargo affinities are uniform and relatively strong ε RR = ε GG = ε RG = 1.7, leading to rapid and simultaneous coalescence of cargo and shell assembly, resulting in a uniform mixture of both cargo types within the complete shell (Multimedia view). In Fig. 3(b), the R-R affinities (ε RR = 1.7) are stronger than the G-G and R-G affinities (ε GG = ε RG = 1.3), leading to simultaneous assembly and coalescence of a nearly pure R domain (Multimedia view).
Figures 3(c)-3(e) show trajectories with a weaker shell-shell interaction of ε SS = 2.0, for which shell nucleation is slow and we observe two-step assembly pathways. In Fig. 3(c), the cargo-cargo affinities are strongest for R-R, and the assembly pathway closely resembles that for a single species-a domain of nearly pure of Chemical Physics shell-cargo, ε SC ; R-R and G-G, ε RR and ε GG ; and R-G, ε RG ). (a) All affinities are relatively strong, driving rapid and simultaneous coalescence of both cargo types and shell assembly (one-step assembly). The encapsulated cargo is uniformly mixed (Multimedia view). (b) At ε SS = 3.5k B T, high ε RR but moderate ε GG and ε RG lead to one-step assembly and coalescence of R cargo, with G cargo essentially excluded from the shell (Multimedia view). (c) At ε SS = 2k B T, strong ε RR but moderate ε GG and ε RG lead to two-step assembly around an almost pure R domain (Multimedia view). (d) Strong ε GG and ε RR and moderate ε RG drive two-step assembly, with coupling between shell closure and cargo compositional fluctuations, leading to encapsulation of nearly pure R and G cargo domains in separate shells. The initial globule that becomes separated is shown in the top cutaway view, while the assembled shells are shown in the bottom view (Multimedia view). (e) Very strong ε RR , strong ε GG , and moderate ε RG lead to two-step assembly around an almost pure G domain. The strength of the ε RR interaction prevents the subunits from successfully encapsulating the globule of R cargo, while the G cargo is able to bud off similarly to (d) and form a properly assembled shell. The cutaway shows the state of the assembled shell and G cargo domain within as well as the R cargo globule that is unable to properly close (Multimedia view R particles condenses, with subsequent assembly and encapsulation by the shell (Multimedia view). We observe a more interesting coupling between cargo composition and shell assembly for strong R-R and G-G affinities (ε RR = ε GG = 1.7) and moderate R-G affinities ε RG = 1.3 [ Fig. 3(d)]. The cargo initially coalesces into a globule that is mixed, albeit with significant compositional fluctuations, consistent with the bulk cargo phase behavior under these parameters [the purple filled hexagon in Fig. 2(b)]. At these parameters, the domain is large enough to form two shells and shell closure couples to compositional fluctuations and occurs preferentially in the vicinity of an interface between R-rich and G-rich domains. In particular, we observe closure at the interface and consequently purified cargo encapsulation, meaning that the encapsulated cargo within an individual shell is nearly pure R or G (Multimedia view).
The origins of coupling between shell closure and compositional fluctuations can be understood from the competing forces that govern shell closure. For the shell to close, curvature of the assembling shell must first generate the formation of a "neck" of cargo, i.e., a narrow region of cargo particles that connects a shell assembly to the remainder of the globule or to other assemblies [e.g., see . The processes of forming and then breaking the neck increase the total cargo interfacial area and correspondingly reduce the number of cargo-cargo interactions. Thus, they are accompanied by a free energy barrier with a height that increases with cargo surface tension and correspondingly increases with cargo-cargo affinity. In Fig. 3(e), the R-R affinities are stronger (ε RR = 2.0), resulting in a larger barrier for closure within the R-rich portion of the globule, and thus a complete shell forms only around G particles (Multimedia view).

Shell closure locks the encapsulated cargo composition out of equilibrium
After a shell closes, fluctuations of subunits large enough to allow cargo particles to escape are extremely rare because all subunits have their maximum number of interactions. Thus, the cargo composition is essentially fixed once the shell completes. However, as noted in Sec. II D, the cargo remains fluid in the range of cargo-cargo affinities that we focus on, and, thus, the cargo can spatially reorganize inside of the shell even after its completion.

D. Factors that control the composition of encapsulated cargo 1. Cargo-cargo affinity
Under most conditions, the cargo-cargo affinities have the strongest effect on the amount and composition of encapsulated cargo. We begin by considering one-step assembly pathways obtained with moderate shell-cargo affinities ε SC = 6.0 and the higher shell-shell affinities ε SS = 3.5. These conditions allow assembly over a broad range of cargo-cargo affinities. most results in this figure with varying R-R and G-G affinities at fixed ε RG = 1.3. First, Fig. 4(c) shows that the amount of packaged cargo increases monotonically with the net encapsulation driving force; i.e., increasing either or both of ε RR and ε GG raises the number of encapsulated cargo particles N C , which asymptotically approaches the value corresponding to crystalline density (N C ≈ 150). 60 This demonstrates that cargo-cargo cohesive interactions are essential to of Chemical Physics obtain full shells. 60 Importantly though, N C does not reach crystalline density until cargo-cargo affinities nearly reach the bulk crystallization threshold ε fs ≈ 3k B T. Reference 60 found that wellformed shells cannot assemble when the cargo crystallizes, because the solid cargo cannot be deformed by assembly subunits, thus preventing shell closure. They also measured radial distribution functions for a single cargo species system, which showed that confinement by the shell imposes layering on the encapsulated cargo due to steric effects, but that the cargo remains fluid. The same effect is present for the bicomponent system studied in the present work, since the two cargo species differ only in their cargo-cargo affinities. Next, to determine the relative loading of the two cargo species, Fig. 4(a) shows the fraction of R species averaged over all shells, f R , FIG. 6. Fraction of shell subunits in complete shells in the two-step assembly pathway regime, ε SS = 2.0 and ε RG = 1.3.

FIG. 7.
Effects of shell-cargo and shell-shell affinity on cargo composition reveal thermodynamic and kinetic influences on encapsulation. (a) The fraction of encapsulated R particles f R as a function of shell-cargo affinity ε SC , for shell-shell affinities ε SS = 2.0, and ε SS = 3.5. To assess the importance of dynamics on these results, f R is also shown for the "equilibrium" simulations in which cargo can exchange between shell interiors and the bulk. Other parameters are ε GG = 1.7 and ε RR = ε RG = 1.3. (b) Cargo composition f R as a function of shell-shell affinity ε SS , at a constant shell-cargo affinity of ε SC = 6.0. Other parameters are ε RR = 1.3, with ε GG = 1.7 and ε RG = 1.3. The "equilibrium" results are also shown, as a straight line since they do not depend on ε SS . as a function of ε GG and ε RR . Further, to characterize the arrangement of cargo particles within shells, Fig. 4(b) shows the fraction of unlike particles in the first solvation shell ( f UL ). The snapshots surrounding Fig. 4(a) illustrate typical arrangements of the encapsulated cargo in each regime. Comparing Figs. 4(a) and 4(b) with results from the bulk cargo system in the absence of microcompartment shells [Figs. 2(b) and 2(e)] shows that the composition and mixing of the two cargo species within shells is similar to the case of bulk cargo coalescence for this parameter regime. However, we will see below that shell assembly has a more significant effect on cargo properties in other parameter regimes. Figure 5 compares these results to two-step assembly pathways that occur for lower shell-shell affinities, ε SS = 2.0. Since we observe significant assembly only for a narrow range of strong cargo-cargo affinities (Fig. 6, discussed next), we present these results as a function of ε RR for fixed ε GG = 1.7. We see that average cargo loading and composition are qualitatively similar between the two cases, although there are quantitative differences revealing deviations from equilibrium, which we consider in more detail in Sec. II D 2.
However, the composition of cargo within individual shells can differ significantly from the mean. Figure 4(e) shows the shellto-shell variability of cargo composition, σ R,shell , defined as the standard deviation of the fraction of R particles within each shell: Because compositional fluctuations depend strongly on the ratio of ε RR and ε GG to ε RG , we present the results as a function of equal R-R and G-G interactions (ε RR = ε GG ), for weak, moderate, and strong ε RG . In each case, as the R-R and G-G affinities exceed ε RG , there is a transition from shells containing uniformly mixed cargo to nearly pure domains of R or G particles respectively, within different shells, corresponding to large σ R,shell .
Interestingly, this variability arises in different ways depending on the class of assembly pathways. To illustrate this difference, Figs. 5(c) and 5(d) show histograms of cargo compositions measured within individual shells for one-step and two-step pathways, respectively. These parameters correspond to the maximum in Fig. 5(a), where f R and f UL are approximately equal for the two pathways, but the mechanisms underlying the cargo separation differ. In one-step The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp pathways, as shown in Fig. 3(b), a compositional fluctuation in the initial cargo globule nucleus leads to preferential coalescence of the same cargo type as shell assembly proceeds, with the shell eventually closing around a nearly pure domain. However, for these parameters, we do occasionally (∼5%) observe shells that encapsulate two domains separated by an interface [see Fig. 5(c)]. In contrast, two-step assembly pathways form purified shells via the compositional fluctuation mechanism discussed in Sec. II C and shown in Fig. 3(d). This leads to highly purified cargo globules within individual shells [see Fig. 5(d)].
a. Cargo-cargo affinities also affect shell assembly. While our focus in this article is on factors that control cargo encapsulation, it is important to note that the driving force for cargo cohesion in turn affects the extent and robustness of shell assembly. Figures 4(d) and 6 show the fraction of subunits in complete shells fs as a function of cargo-cargo affinities for ε SS = 2.0 and ε SS = 3.5, respectively. We observe that assembly occurs over a broad range of cargo-cargo affinity values for the higher shell-shell affinity [ε SS = 3.5, Fig. 4(d)], but for the lower shell-shell affinity (ε SS = 2.0, Fig. 6), significant assembly requires at least one strong cargo-cargo affinity; i.e., εij ≥ 1.7 for at least one of ij ∈ {RR, GG, RG}. This result indicates that nucleation of shell assembly for ε SS ≲ 2.5 requires adsorption onto a cargo globule and thus requires that cargo-cargo affinities are sufficient to drive cargo phase separation. This observation is consistent with results from simulations with a single cargo species. 60 2. Shell-shell and shell-cargo affinities affect cargo packaging through equilibrium and nonequilibrium mechanisms Although the cargo-cargo affinities most strongly affect the amount and properties of encapsulated cargo, the shell protein affinities (ε SC and ε SS ) can also modulate cargo encapsulation. These effects depend on an interplay between thermodynamic and kinetic effects. The assembly trajectories in our simulations begin from an initial condition of dispersed subunits and cargo that are significantly out of equilibrium. Moreover, once a shell closes, the timescale for reconfiguration to another shell morphology or exchange of the encapsulated cargo is much longer than assembly timescales since all subunits have their maximum number of interactions. Therefore, even though the dynamics satisfy microscopic reversibility, trajectories are not guaranteed to approach equilibrium on the finite timescales of simulations and experiments. 157 To distinguish between equilibrium and kinetic effects, we compare the results of our assembly simulations to the alternative "equilibrium" or permeable-shell system described above, in which shells are completely assembled but are made permeable to cargo particles, thus allowing the encapsulated cargo to fully equilibrate with cargo particles in solution (see Sec. IV). The amount and composition of encapsulated cargo in these simulations is, thus, independent of assembly kinetics (and shell-shell affinities). Figure 8 compares the results of assembly trajectories at parameters leading to one-step assembly (corresponding to Fig. 4) with the equilibrium simulations. Over much of parameter space, the amount and compositions of encapsulated cargo are similar between the two systems, indicating that trajectories are near equilibrium during assembly. This can be understood since the diffusional exchange of cargo particles between the partially encapsulated cargo globule and bulk is rapid compared to shell growth under those conditions. However, we do observe significant deviations from equilibrium cargo encapsulation under certain parameter ranges. These effects are strongly dependent on the shell-shell and shell-cargo affinities, since they are key determinants affecting the extent to which assembly is out of equilibrium.
For example, Fig. 5(a) shows the dependence of f R on ε RR for one-step (ε SS = 3.5) and two-step (ε SS = 2.0) assembly pathways. For moderate cargo-cargo affinities ε RR ≲ 1.5, the two-step assembly and equilibrium results closely agree, while the one-step simulations result in higher f R . This difference arises because the initial globule is enriched in R particles (since the globule nucleation is driven in part by shell-cargo attractive interactions that are agnostic to particle type), and the rapid one-step assembly pathways do not allow the cargo globule time to fully equilibrate with the bulk. In contrast, the slow onset of nucleation and assembly in the two-step pathways allows globule equilibration. Interestingly, the situation reverses at high ε RR ≳ 1.6, where there is an abrupt drop in f R for the two-step case. This arises because the strong R-R affinities favor segregation of the two cargo species and a high surface tension that inhibits closure of a shell around R particles. A comparable trajectory is shown in Fig. 3(e). Notice that a shell successfully assembles around the cargo particles with weaker affinities (G), but only a malformed structure forms around the R globule because the shell subunits are unable to close around a commensurate-sized piece of the globule. Note that we do not show the permeable-shell results for such high ε RR because the strong cargo-cargo affinities result in a cargo globule that extends outside of the shell, which is an artifact of the permeable-shell condition imposed in those simulations.
We also observe that the extent of cargo mixing within shells [ Fig. 5(b)] decreases with increasing ε RR . This is partially an The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp equilibrium effect due to the increasing drive for cargo segregation, especially as the R-R and G-G affinities exceed the R-G affinity (ε RG = 1.3). However, while the two-step results are close to equilibrium, mixing is further reduced for the one-step assembly trajectory. This effect can be attributed to the coupling between composition fluctuations and interfacial effects within the globule during the initial stage of the two-step pathway described in Sec. II C and Fig. 3(d). Due to interfacial tension, the boundary between the two domains provides a "weak point" for insertion of shell subunits, and, thus, shell closure tends to occur in the vicinity of the interface. The closing shell becomes enriched in the cargo species that dominates the domain on the same side of the interface as the assembling shell; this composition is then effectively locked in by shell closure.
Having compared the differences in cargo encapsulation between the two classes of assembly pathways, we now consider parameter values between these two extremes. Figs. 7(a) and 7(b), respectively, show the effects of shell-cargo and shell-shell affinities on the composition of encapsulated cargo. Here, we have focused on the case for which cohesive interactions between pairs of G particles are stronger than for R-R and R-G, ε GG = 1.7 and ε RR = ε RG = 1.3, for which the equilibrium and bulk composition of R particles is about 20%. Fig. 7(a) shows that for weak shell-cargo affinities, the encapsulated fraction of the R particles is significantly smaller than the equilibrium value, but that it monotonically increases with ε SC . The trend is the same for both shell-shell affinities, except that higher shell-cargo affinities are required to observe assembly for lower shell-shell affinities. Similarly, Fig. 7(b) shows that f R increases monotonically with shell-shell affinity. These trends can be understood as follows. At low ε SC , the cargo-cargo affinities are the primary driving force for encapsulation, which at these parameters favors encapsulation of the G cargo and, thus, small f R . The fact that f R in the dynamical assembly simulations is smaller than that observed in the bulk cargo and equilibrium simulations suggests that the initial globule that nucleates is enriched in G, and that the shell closes before there is sufficient time for the globule composition to equilibrate with the bulk. However, as ε SC or ε SS increases, the strong shell-cargo affinities enhance nucleation, allowing a cargo globule to form with a composition closer to the bulk composition, f R = 0.5.

III. CONCLUSIONS
We have performed equilibrium and dynamical simulations to investigate the encapsulation of multicomponent fluid cargoes by self-assembling microcompartments. To elucidate mechanisms that control the amount and composition of cargo encapsulation, we compared these results against simulations of cargo coalescence in the absence of microcompartment assembly.
Our simulation results show that the binding affinities among the different cargo species are the strongest determinant of the assembly pathway as well as the amount, composition, and spatial organization of cargo within assembled shells. Analogous to the case of encapsulation of a single component cargo, 57,60,61,120,122,123,161 the simulations exhibit two broad classes of assembly pathways. Relatively strong cargo-cargo affinities and low shell-shell affinities lead to two-step assembly in which the cargo first forms a globule, followed by adsorption and assembly of the shell around the globule. Weaker cargo-cargo affinities and higher shell-shell affinities drive one-step assembly in which cargo coalescence and shell assembly occur simultaneously. However, the bicomponent cargo allows for more diverse assembly pathways and morphologies, depending on when and if the two components mix within the encapsulated globule. In particular, we observe assembly pathways ranging from coalescence and encapsulation of a nearly pure globule of a single cargo species, to coalescence of a phase-separated globule with each phase becoming encapsulated in different shells, to encapsulation of a globule consisting of a uniform mixture of both species. Correspondingly, the assembled shells can contain a variety of cargo compositions and morphologies, ranging from highly enriched in a single species to nearly uniform mixing of the two species. In regimes where both cargo species are encapsulated, the two species may be uniformly mixed within each shell, or separated, with each individual shell containing an encapsulated cargo that is nearly pure in one or the other species.
While the relative strengths of the different cargo-cargo binding affinities (ε RR , ε GG , and ε RG ) determine cargo composition and mixing, the net cargo-cargo affinity most strongly influences the amount of cargo that is encapsulated. In particular, for low cargo-cargo affinities, we observe poor cargo loading (e.g., shells that are nearly empty except for a layer of cargo at the surface). Further, depending on the strength of the subunit-subunit affinities, there is a threshold value of cargo-cargo affinities below which assembled microcompartments are either unstable or fail to nucleate on relevant timescales [Figs. 4(d) and 6].
While cargo composition and mixing within assembled shells is qualitatively similar to the behavior of cargo coalescence in bulk (in the absence of microcompartments) in many parameter regimes, our simulations also identify parameter regimes in which the shell-cargo and shell-shell binding affinities can be manipulated to control cargo encapsulation (e.g., Figs. 5 and 7). A combination of thermodynamic and kinetic mechanisms underlie these effects. In particular, we see a marked difference between the properties of the encapsulated cargo in shells that assemble by one-step or twostep pathways. When cargo-cargo affinities are asymmetric, leading to preferential encapsulation of one species, weaker shell-cargo and shell-shell affinities and two-step assembly pathways tend to enhance the degree of preferential encapsulation, leading to more purified cargo. This behavior reflects the fact that the cargo globule that nucleates in the first step tends to be enriched in the species with higher affinities, and, if assembly occurs more rapidly than the process of cargo exchange with bulk, shell closure locks in this outof-equilibrium composition. Similarly, when the set of cargo-cargo affinities places the system close to ternary phase coexistence (a dilute phase and two dense phases, respectively, rich in either cargo species), two-step assembly pathways can enhance separation of the cargo species. This result can be explained at least in part by the tendency of assembling shells to close in the vicinity of transient interfaces between domains of the two cargo species (Fig. 2).
The latter predictions can be tested in microcompartment systems by performing mutagenesis of shell protein-protein binding interfaces (to alter shell-shell affinities) or "encapsulation peptides" that mediate shell-cargo interactions. [162][163][164][165] Similarly, cargo-cargo affinities may be tuned through mutagenesis, but they are less well understood at present. Depending on the system, cargo-cargo affinities may arise from direct pair interactions, scaffold-mediated interactions, or a combination of the two. 31,36,39,121,122,126,166,167 The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp

A. Outlook
It will be interesting to extend our model to explicitly incorporate scaffold-mediated cargo-cargo and shell-cargo interactions, for example as Mohajerani et al. 62 did for a single-component cargo. More broadly, the predictions from our models can be extended beyond microcompartments, for example to reengineered viral capsids [46][47][48]51,[147][148][149][150][151][152][153][154][155] or synthetic capsids constructed via protein engineering (e.g., Refs. [168][169][170] or DNA origami 171 to assemble around specific cargoes. In some of these cases, the shell-cargo interactions are driven by "scaffolds" attached to the inner surface of the capsid. For example, many virus capsid proteins have flexible peptide tails that drive RNA encapsulation, 172-174 while Sigl et al. 171 designed DNA origami capsid subunits with single-stranded DNA oligomers that bound to complementary strands on cargo particles. Previous simulations of microcompartments 62 and virus assembly 106,[175][176][177][178] suggest that the properties of such scaffolds can be designed to manipulate the spatial organization of the cargo. Finally, it will be important to allow for shell assembly into different geometries to study coupling between shell size and cargo composition.

A. Computational model
We consider a minimal model for microcompartment assembly and encapsulation of two species of cargo particles, denoted as "R" and "G." Microcompartment shells assemble from pentameric, hexameric, and pseudo-hexameric (trimeric) protein oligomers [e.g., Fig. 3(a) in Ref. 58 and Refs. 18,19,and 179]. Experiments suggest that these oligomers are the basic assembly units, meaning that smaller protein complexes do not contribute significantly to the assembly process. 19,180 To focus on factors that control the composition of encapsulated cargo, we consider a minimal shell model developed in Perlmutter et al., 60 with only pentamer and hexamer subunits that have interactions designed to drive assembly into a T = 3 shells (containing 12 pentamers and 20 hexamers in a truncated icosahedron geometry). Our model builds on previous simulations of microcompartment assembly with no cargo 66 or one cargo species, 60,61,65 as well as previous models for virus assembly. 76,79,80,83,102,[175][176][177][181][182][183]

Shell-shell interactions
In our model, shell subunits interact through both repulsive and attractive forces. The repulsions consist of "Top" (T) pseudoatoms that exist above the plane of the hexamer and "Bottom" (B) pseudoatoms that exist below the plane. More specifically, we denote the top pseudoatoms as "TH" and "TP," respectively, for hexamer and pentamer subunits, and similarly the bottom pseudoatoms as "BH" and "BP." All pairs of Top and Bottom pseudoatoms (excepting those on the same subunit) interact via a repulsive Lennard-Jones potential [Eq. (B3)]. The Top-Top diameters and repulsion strengths set the shell spontaneous curvature and bending modulus (Appendix B).
Microcompartment protein shell assembly is primarily driven by interactions along the edges of the hexameric and pentameric subunits. 58 We represent these interactions by "Attractors" on the perimeter of each shell subunit. Complementary Attractors on nearby subunits have short-range interactions modeled by a Morse potential [Eq. (B4) in Appendix B]. Attractors that are not complementary do not interact. The arrangement of Attractors on subunit edges is shown in Fig. 1, with pairs of complementary Attractors indicated by cyan double-headed arrows. The shell-shell binding affinity is proportional to the well-depth of the Morse potential between complementary Attractors, ε SS .
To focus on the role of interaction strengths in determining cargo encapsulation, we consider shell subunits that preferentially form only T = 3 shells. To enforce this restriction, only three of the six edges of a hexamer have attractive interactions with pentamers, and there are no pentamer-pentamer attractive interactions (see Fig. 1).

Shell-cargo interactions
We model attractive interactions between hexamer subunits and the cargo particles by a Morse potential between cargo particles and subunit "B" pseudoatoms, with well-depth ε SC . These interactions represent the effect of "encapsulation peptides" that target cargo to microcompartment interiors by mediating cargo-hexamer interactions. 31,39,121,122,166 To minimize the number of parameters, we keep ε SC the same for both cargo species. We also add a layer of "Excluders" in the plane of the "Top" pseudoatoms, which represent shell-cargo excluded volume interactions. Since the shell-shell interaction geometries are already controlled by the Attractor, Top, and Bottom pseudoatoms, we do not consider Excluder-Excluder interactions.

Cargo-cargo interactions
In natural microcompartment systems, the interior cargo undergoes phase separation prior to or during shell assembly. The attractions between cargo particles that drive phase separation may be mediated by microcompartment scaffolding proteins (e.g., Refs. 120, 122, and 126) and/or direct pair interactions between cargo particles. 36 Similarly, in synthetic systems, it is possible to engineer direct pair or scaffold-mediated attractions cargo molecules. In our model, we represent these scenarios in a minimal manner. We model both species of cargo particles with spherically symmetric excluded volume and direct pair attractions, implemented via an attractive Lennard-Jones (LJ) potential, with welldepth values ε GG , ε RR , and ε RG for pairs of G-G, R-R, and R-G particles.

B. Mapping simulation parameters to biological microcompartments 1. Relevant ranges of affinities
Although affinities among microcompartment components have not been measured, we estimate relevant ranges of our control parameters based on experimental observations of carboxysomes as follows: a. Shell-shell affinities. Because empty carboxysome shells are not typically observed in cells (for wild-type carboxysomes), we have focused on shell-shell affinities that are below the threshold for empty shell assembly (ε SS ≈ 4k B T). Furthermore, experimental observations suggest that α-carboxysomes assemble by one-step assembly pathways, 57,120,161 whereas β-carboxysomes 122,123 and many other microcompartments 23,179 assemble by two-step The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp pathways. Therefore, except where mentioned otherwise, we consider two cases (ε SS /k B T = 2, 3.5), which are, respectively, below and above the transition between one-step and two-step assembly pathways (discussed in Sec. II C).
b. Shell-cargo affinities. Since the shell subunits robustly assemble around RuBisCO in both α and β carboxysomes, we have focused on shell-cargo affinities that are strong enough to observe this (ε SC = 6k B T, free energy g SC ≈ −14k B T) except where mentioned otherwise.
c. Cargo-cargo affinities. For the case of α-carboxysomes, the RuBisCO does not undergo phase separation in the absence of shell proteins at physiological conditions, 57,120,161 but it can phase separate without shell proteins under nonphysiological conditions. 126,128 Thus, we infer that the effective cargo-cargo affinity for α-carboxysomes is near but below the threshold for phase separation, which corresponds to a cargo-cargo affinity of ε vl ≈ 1.3k B T at our simulated cargo concentration. In contrast, for β-carboxysomes, the RuBisCO undergoes phase separation before shell assembly, suggesting a cargo-cargo affinity above ε vl but below crystallization (ε fs ). Our simulations span this range.

Length scales
Although the model is designed to be generic, we can approximately map it to carboxysomes by setting the cargo diameter (the unit length scale in the model) to that of the RuBisCO holoenzyme, implying r * ≈ 13 nm. However, to enable tractable simulation of long assembly timescales, we have set the subunit/cargo size ratio to be larger than in carboxysomes-model subunits have a side length of r * and are thus about three times larger than carboxysome hexamers (side length ≈4 nm). See Ref. 62 for further discussion.

C. Simulations and systems
We performed simulations using the Langevin dynamics algorithm in HOOMD [which uses graphical processing units (GPUs) to accelerate molecular dynamics simulations 184 ], and periodic boundary conditions in a cubic box to represent a bulk system. The subunits were modeled as rigid bodies. 185 Each simulation was performed in the NVT ensemble, using the HOOMD fundamental units, 186 with the unit length scale r * defined as the circumradius of the pentagonal subunit (the cargo diameter is also set to r * ), and energies given in units of the thermal energy, k B T. The simulation timestep was 0.005 in dimensionless time units for all systems.

Systems
We simulated several systems as follows: a. Shell assembly and cargo encapsulation. For dynamical simulations of shell assembly in the presence of cargo, each simulation contained enough subunits to form four complete microcompartments (48 pentamers and 80 hexamers), in a cubic simulation box with side length 40r * . Each simulation also contained 854 cargo particles, with composition 50% "R" and "G," respectively. With r * = 13 nm, this corresponds to concentrations of cp = 0.6, c h = 1, and cc = 10 μM for pentamers, hexamers, and cargo.
We performed most simulations for 2 × 10 6 timesteps; simulations in the two-step assembly regime were performed for longer, 2 × 10 7 timesteps, as they assemble more slowly.
b. Bulk simulations without shell subunits. We also performed large simulations of cargo without shell subunits, to determine the phase behavior of the cargo as a function of the cargo-cargo affinities. Each of these simulations contained 46 938 cargo particles (50% "R" and "G," respectively) in a box with side length 160, to maintain the same cargo concentration as in the shell assembly simulations. We performed each simulation for 2 × 10 7 timesteps. c. Equilibrium cargo encapsulation. To identify effects of outof-equilibrium dynamics of shell assembly on the properties of the encapsulated cargo, we performed additional simulations in which we simulated cargo dynamics in the presence of an assembled shell that was made permeable to cargo. In particular, these simulations consisted of a fully assembled immobile shell placed in the center of the box, with one missing pentamer to allow permeability to occur with the surrounding cargo. While the shell could not move, the ε SC was still active. Each simulation contained 854 cargo particles with box side length 40r * to maintain the same cargo concentration and composition as in the other systems. We performed each simulation for 2 × 10 6 timesteps, which we determined was sufficient to achieve equilibration of cargo inside of the shell by monitoring the amount, composition, and mixing of the encapsulated cargo.

Initial conditions
Except for the equilibrium cargo encapsulation systems, all simulations were initialized with random positions for cargo particles and shell subunits (except that significant overlaps between pseudoatoms were forbidden). The equilibrium cargo encapsulation simulations were initialized with an assembled microcompartment shell whose subunits were fixed throughout the simulation, and the cargo particles were initialized with random locations with overlaps forbidden.

Statistical significance
For most systems, we performed ten independent trials at each parameter set. For most parameter sets at which we report statistical errors, we performed 20 independent trials. For parameter sets that resulted in low probability events (≤5%), we performed 60 trials to ensure low statistical error.

Conflict of Interest
The authors have no conflicts to disclose.

DATA AVAILABILITY
The data that support the findings of this study are openly available in Open Science Framework at http://doi.org/10.17605/ OSF.IO/SKZ7F.

APPENDIX A: SHELL-SHELL AND SHELL-CARGO BINDING FREE ENERGY VALUES
We estimated the subunit-subunit binding free energy values as functions of the well-depth parameter ε SS by measuring the dimerization equilibrium constant in simulations of subunits only capable of forming dimers. For both pentamer-hexamer and hexamer-hexamer dimers, we obtain binding free energies that are linear functions of the well-depth ε SS . 60 We interpret the y-intercept as the binding entropy. Setting r * = 13 nm and using the conventional standard state concentration of 1M (rather than r 3 * as in Ref. 60) gives where g ph and g hh are the binding affinities for pentamer-hexamer and hexamer-hexamer pairs. The dissociation constants are then given by K hh d = exp(g hh /k B T) and K ph d = exp(g ph /k B T). Table I shows the results for hexamer-hexamer dimerization; the values for pentamer-hexamer association are approximately the same. Similarly, we estimated the shell subunit adsorption free energy by performing simulations of subunits that cannot assemble (ε SS = 0) in the presence of a cargo globule as 60 For the free energy of shell assembly, we consider a shell comprised of npent = 12 pentamers and n hex hexamers, which have n ph pentamer-hexamer contacts with binding energy ε ph and n hh hexamer-hexamer contacts with energy ε hh . For our T = 3 model, n hex = 20, n ph = 60, and n hh = 30. The assembly free energy is then given by ΔG shell = n ph ε ph + n hh ε hh − T(npentspent + n hex s hex + s config ) where s config accounts for the configurational entropy associated with subunit and shell symmetries: s config = k B ln(5 n pent 3 n hex /60), and γ and A C are the surface tension and area of the cargo globule. Finally, we define an effective dissociation constant for the complex as K shell d ≡ k B T ln Δg shell with Δg shell = ΔG shell /(npent + n hex ) the free energy per subunit averaged over pentamers and hexamers (see Table I). However, this calculation is only approximate since it neglects changes in binding entropy between dimerization and the full complex. Furthermore, to simplify the calculation and focus on the shell subunits, we have neglected cargo-cargo contributions to the complex stability. This provides a reasonable estimate for cargo-cargo affinities near liquid-vapor coexistence, but it over/underestimates the shell stability for weaker/stronger cargo-cargo affinities. For a full analysis including cargo-cargo contributions, see Perlmutter et al. 60

APPENDIX B: MODEL DETAILS
In our model, all potentials can be decomposed into pairwise interactions. Potentials involving shell subunits further decompose into pairwise interactions between their constituent building blocks-the Excluders, Attractors, "Top," and "Bottom" pseudoatoms. It is convenient to express the total energy of the system as the sum of three terms, involving shell-shell (U SS ), cargo-cargo (U CC ), and shell-cargo (U SC ) interactions, each summed over all pairs of the appropriate type, where ∑ shell i ∑ shell j<i is the sum over all distinct pairs of shell subunits in the system, ∑ shell i ∑ cargo j is the sum over all shell-cargo particle pairs, etc.

Shell-shell interaction potentials
The shell-shell potential U SS is the sum of the attractive interactions between complementary Attractors, and geometry guiding repulsive interactions between "Top"-"Top," "Bottom"-"Bottom," and "Top"-"Bottom" pairs. There are no interactions between members of the same rigid body. Thus, for notational clarity, we index rigid bodies and nonrigid pseudoatoms in Roman, while the pseudoatoms comprising a particular rigid body are indexed in Greek. For subunit i, we denote its Attractor positions as {aiα} with the set comprising all Attractors α, its "Top" position ti, and its "Bottom" position bi.
The shell-shell interaction potential between two subunits i and j is then defined as U SS ({aiα}, ti, aj, tj) = ε SS Ł(|ti − tj|, σt,ij) + ε SS Ł(|bi − bj|, σ b ) + ε SS Ł(|bi − tj|, σ tb ) + N ai ,N aj ∑ α,β ε SS ℳ (|a iα − a jβ |, r 0 , ρ, r att cut ). (B2) The function Ł is defined as the repulsive component of the Lennard-Jones potential shifted to zero at the interaction diameter, with V shift (rcut) being the value of the (unshifted) potential at rcut. The parameter ε SS sets the strength of the shell-shell attraction at each Attractor site, Nai is the number of Attractor pseudoatoms in subunit i, and ε angle scales the repulsive interactions that enforce the geometry.

Shell-shell interaction parameter values a. Attractors
The strength of attractive interactions is parameterized by the well-depth ε SS for a pair of Attractors on hexamers as follows. Hexamer-hexamer edge Attractor pairs (A2-A6, A3-A5, and A5-A6) have a well-depth of ε SS . Because vertex Attractors (A1, A4) have multiple partners in an assembled structure, whereas edge Attractors have only one, the well-depth for the vertex pairs (A1-A4 and A4-A4) is set to 0.5ε SS . Similarly, for pentamer-hexamer interactions, the well-depth for edge Attractor pairs (A2-A5, A3-A6) is ε SS , while the vertex interaction pairs (A1-A4 and A4-A4) have 0.5ε SS .

b. Repulsive interactions
The "Top" and "Bottom" heights, or distance out of the Attractor plane, are set to h = 1/2r b , with r b = 1 the distance between a vertex Attractor and the center of the pentagon. For all simulations, shells have T = 3 preferred curvature and σ tb = 1.8r b is the diameter of the "Top"-"Bottom" interaction (this prevents subunits from binding in inverted configurations 76 ) while σ b = 1.5r b is the diameter of the "Bottom"-"Bottom" interaction. In contrast to the latter parameters, σt,ij, the effective diameter of the "Top"-"Top" interaction, depends on the species of subunits i and j; denoting a pentagonal or hexagonal subunit as "p" or "h," respectively, σt, pp = 2.1r b , σ t, hh = 2.4r b , and σ t, ph = 2.2r b . The parameter r 0 is the minimum energy Attractor distance, set to 0.2r b , ρ = 4r b determines the width of the attractive interaction, and r att cut = 2.0r b is the cutoff distance for the Attractor potential. Since the interactions just described are sufficient to describe assembly of the shell subunits, we included no Excluder-Excluder interactions.

Cargo-cargo interactions
The interaction between cargo particles is given by a sum over all interacting pairs, with ℒ being the full Lennard-Jones interaction, and N RR , N GG , and N RG as the number of R-R, G-G, and R-G pairs in the system, and the cargo diameter σ C = r b and cutoff r C cut = 3σ C are the same for all cargo species.

Shell-cargo interactions
The shell-cargo interaction is modeled by a short-range repulsion between cargo-Excluder and cargo-"Top" pairs representing the excluded volume, plus an attractive interaction between pairs of cargo particles and shell subunit "Bottom" pseudoatoms. For subunit i with Excluder positions {xiα} and "Bottom" psuedoatom bi, and cargo particle j with position Rj, the potential is given by where ε SC parameterizes the shell-cargo interaction strength, Nx, Nt, and N b are the numbers of Excluders and "Top" and "Bottom" pseudoatoms on a shell subunit, σex = 0.5r b and σt = 0.5r b are the effective diameters of the Excluder-cargo and "Top"-cargo repulsions, r SC 0 = 0.5r b is the minimum energy Attractor distance, the width parameter is ρ SC = 2.5r b , and the cutoff is set to r SC cut = 3.0r b .